%
% @brief 射频系统的 冲激响应特性 补偿
% @author 41204@qq.com / 2025-10
% @note 本类中的函数均以 Static 方式定义，以方便使用。
%

classdef rf_hfreq_compensate
    methods(Static)
        %% 补偿射频接收机的幅度和相位失真
        % 输入：
        %   A: 输入IQ信号 (复数)
        %   B: 输出IQ信号 (复数)  
        %   fs: 采样频率
        % 输出：
        %   A_compensated: 补偿后的信号
        %   H_estimated: 估计的系统频率响应
        %   compensation_results: 补偿效果分析
        function [A_compensated, H_estimated, compensation_results] = compensate_rf_response(A, B, fs, varargin)        
            p = inputParser;
            addParameter(p, 'smoothing_factor', 0.01, @(x) x>=0 && x<=1);
            addParameter(p, 'compensation_limit_db', 20, @isnumeric);
            parse(p, varargin{:});
            
            N = length(A);
            
            % 1. 估计系统频率响应
            A_fft = fft(A, N);
            B_fft = fft(B, N);
            
            % 计算原始频率响应
            H_raw = B_fft ./ A_fft;
            
            % 处理除零和噪声问题
            small_values = abs(A_fft) < max(abs(A_fft)) * 1e-6;
            H_raw(small_values) = 1;  % 在这些频率点设为无失真
            
            % 2. 平滑频率响应（重要！）
            H_smoothed = rm.rf_hfreq_compensate.smooth_frequency_response(H_raw, p.Results.smoothing_factor);
            H_estimated = H_smoothed;
            
            % 3. 设计补偿滤波器（逆系统）
            H_comp = 1 ./ H_estimated;
            
            % 限制最大补偿幅度（避免噪声放大）
            max_comp_mag = 10^(p.Results.compensation_limit_db/20);
            H_comp_mag = abs(H_comp);
            H_comp(H_comp_mag > max_comp_mag) = max_comp_mag .* exp(1j * angle(H_comp(H_comp_mag > max_comp_mag)));
            
            % 4. 应用补偿
            B_fft_comp = B_fft .* H_comp;
            A_compensated = ifft(B_fft_comp, 'symmetric');
            
            % 5. 评估补偿效果
            compensation_results = rm.rf_hfreq_compensate.evaluate_compensation(A, A_compensated, B);
            
            % 6. 可视化结果
            rm.rf_hfreq_compensate.plot_compensation_results(A, B, A_compensated, H_estimated, fs, compensation_results);
        end
        
        function H_smoothed = smooth_frequency_response(H_raw, smoothing_factor)
            % 使用移动平均或低通滤波平滑频率响应
            if smoothing_factor > 0
                window_size = max(3, round(length(H_raw) * smoothing_factor));
                if mod(window_size, 2) == 0
                    window_size = window_size + 1;
                end
                
                % 幅度平滑
                mag_smoothed = movmean(abs(H_raw), window_size);
                
                % 相位平滑（需要小心处理相位卷绕）
                phase_raw = angle(H_raw);
                phase_unwrapped = unwrap(phase_raw);
                phase_smoothed = movmean(phase_unwrapped, window_size);
                
                H_smoothed = mag_smoothed .* exp(1j * phase_smoothed);
            else
                H_smoothed = H_raw;
            end
        end
        
        function results = evaluate_compensation(A_original, A_compensated, B_output)
            % 评估补偿效果
            results.original_mse = mean(abs(A_original - B_output).^2);
            results.compensated_mse = mean(abs(A_original - A_compensated).^2);
            results.improvement_db = 10 * log10(results.original_mse / results.compensated_mse);
            
            results.original_correlation = corr(real(A_original), real(B_output));
            results.compensated_correlation = corr(real(A_original), real(A_compensated));
            
            results.evml_original = rm.rf_hfreq_compensate.calculate_evm(A_original, B_output);
            results.evml_compensated = rm.rf_hfreq_compensate.calculate_evm(A_original, A_compensated);
        end
        
        function evm = calculate_evm(reference, measured)
            % 计算误差向量幅度
            error = reference - measured;
            evm = sqrt(mean(abs(error).^2) / mean(abs(reference).^2)) * 100;
        end

        %% 示例：测试射频响应补偿
        function test_rf_compensation()
            fs = 100e6;  % 100 MHz 采样率
            t = (0:9999)' / fs;
            
            % 1. 生成测试信号
            % 使用QPSK信号作为测试信号
            symbols = (2*randi([0 1], 1000, 1)-1) + 1j*(2*randi([0 1], 1000, 1)-1);
            A_original = rectpulse(symbols, 10);  % 上采样
            A_original = A_original(1:length(t));  % 截断到合适长度
            
            % 2. 模拟射频接收机失真
            % 幅度失真（带通特性）
            f = linspace(-fs/2, fs/2, length(t));
            magnitude_response = 1 - 0.3 * exp(-(f/20e6).^2);  % 高斯形状的幅度响应
            
            % 相位失真（非线性相位）
            phase_response = 0.1 * (f/10e6).^3;  % 立方相位响应
            
            % 组合频率响应
            H_system = magnitude_response .* exp(1j * phase_response);
            H_system = ifftshift(H_system);  % 转换到FFT顺序
            
            % 3. 应用系统失真
            A_fft = fft(A_original);
            B_fft = A_fft .* H_system.';
            B_output = ifft(B_fft, 'symmetric');
            
            % 添加噪声模拟真实环境
            SNR_dB = 30;
            noise_power = var(B_output) / (10^(SNR_dB/10));
            B_output = B_output + sqrt(noise_power/2) * (randn(size(B_output)) + 1j*randn(size(B_output)));
            
            % 4. 应用补偿算法
            [A_compensated, H_estimated, results] = rm.rf_hfreq_compensate.compensate_rf_response(A_original, B_output, fs, 'smoothing_factor', 0.02, 'compensation_limit_db', 15);
            
            % 5. 显示结果
            fprintf('补偿效果评估:\n');
            fprintf('原始 MSE: %.6f\n', results.original_mse);
            fprintf('补偿后 MSE: %.6f\n', results.compensated_mse);
            fprintf('改善程度: %.2f dB\n', results.improvement_db);
            fprintf('原始 EVM: %.2f%%\n', results.evml_original);
            fprintf('补偿后 EVM: %.2f%%\n', results.evml_compensated);
        end
        
        function plot_compensation_results(A, B, A_comp, H_estimated, fs, results)
            % 绘制补偿结果
            figure('Position', [100, 100, 1200, 800]);
            
            % 时域信号对比
            subplot(2,3,1);
            plot(real(A(1:200)), 'b-', 'LineWidth', 2); hold on;
            plot(real(B(1:200)), 'r--', 'LineWidth', 1);
            plot(real(A_comp(1:200)), 'g:', 'LineWidth', 1.5);
            legend('原始输入', '失真输出', '补偿后');
            title('时域信号对比 (实部)');
            xlabel('样本点'); ylabel('幅度');
            grid on;
            
            % 星座图对比
            subplot(2,3,2);
            plot(real(B), imag(B), 'r.', 'MarkerSize', 1); hold on;
            plot(real(A_comp), imag(A_comp), 'g.', 'MarkerSize', 1);
            plot(real(A), imag(A), 'bo', 'MarkerSize', 2, 'LineWidth', 1.5);
            legend('失真输出', '补偿后', '原始输入');
            title('星座图对比');
            axis equal; grid on;
            
            % 估计的频率响应
            subplot(2,3,3);
            f_axis = (-length(H_estimated)/2:length(H_estimated)/2-1) * fs / length(H_estimated);
            plot(f_axis/1e6, 20*log10(abs(fftshift(H_estimated))));
            title('估计的系统幅度响应');
            xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
            grid on;
            
            % 误差统计
            subplot(2,3,4);
            errors = [results.original_mse, results.compensated_mse;
                      results.evml_original, results.evml_compensated];
            bar(errors');
            set(gca, 'XTickLabel', {'MSE', 'EVM (%)'});
            legend('补偿前', '补偿后');
            title('误差指标对比');
            grid on;
            
            % 相关函数
            subplot(2,3,5);
            [corr_orig, lags] = xcorr(real(A), real(B), 50, 'coeff');
            [corr_comp, ~] = xcorr(real(A), real(A_comp), 50, 'coeff');
            plot(lags, corr_orig, 'r-', lags, corr_comp, 'g-');
            legend('与失真信号相关', '与补偿信号相关');
            title('互相关函数');
            grid on;
            
            % 频谱对比
            subplot(2,3,6);
            [P_orig, f] = pwelch(A, 256, 128, 256, fs);
            [P_dist, ~] = pwelch(B, 256, 128, 256, fs);
            [P_comp, ~] = pwelch(A_comp, 256, 128, 256, fs);
            plot(f/1e6, 10*log10(P_orig), 'b-', f/1e6, 10*log10(P_dist), 'r--', f/1e6, 10*log10(P_comp), 'g:');
            legend('原始输入', '失真输出', '补偿后');
            title('功率谱密度对比');
            xlabel('频率 (MHz)'); ylabel('PSD (dB/Hz)');
            grid on;
        end
    end
end
